{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a7638400",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#吉布斯采样过程\" data-toc-modified-id=\"吉布斯采样过程-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>吉布斯采样过程</a></span></li><li><span><a href=\"#特征筛选\" data-toc-modified-id=\"特征筛选-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>特征筛选</a></span></li><li><span><a href=\"#评价指标对比\" data-toc-modified-id=\"评价指标对比-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>评价指标对比</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4e443b9",
   "metadata": {},
   "source": [
    "# 吉布斯采样过程"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dcb1610",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import random\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import datetime as dt\n",
    "import statsmodels.api as sm\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model import LogisticRegression \n",
    "from sklearn.metrics import cohen_kappa_score\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.metrics import roc_auc_score\n",
    "from sklearn.metrics import roc_curve\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "matplotlib.rcParams['font.sans-serif']=['SimHei']   \n",
    "matplotlib.rcParams['axes.unicode_minus']=False "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff578bfe",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gibbs_sampling(data,sample,turns,temp,criterion):\n",
    "    endog = sample.fillna(0)\n",
    "    def ic(df,criterion):\n",
    "        exog = sm.tools.add_constant(df)\n",
    "        m = sm.Logit(endog.loc[exog.index],exog)\n",
    "        reg = m.fit(maxiter = 50)\n",
    "        if criterion == 'bic':\n",
    "            return reg.bic\n",
    "        elif criterion == 'aic':\n",
    "            return reg.aic\n",
    "\n",
    "    # ------------------------- \n",
    "    record = {}\n",
    "    features = data.columns.tolist()\n",
    "    result = {}\n",
    "    criterion_value = {}\n",
    "    \n",
    "    # -------------------------------------------------------------------------\n",
    "    init = pd.Series(np.random.rand(len(features)), index = features)\n",
    "    init = init.where(init >= 15 / len(features), 1).where(init < 15 / len(features), 0).to_dict()\n",
    "    result[0] = init\n",
    "    for i in range(1, turns + 1):\n",
    "        result[i] = result[i - 1].copy()\n",
    "        record[i] = {}\n",
    "        criterion_value[i] = {}\n",
    "        for j in random.sample(features, len(features)):\n",
    "            k = 1\n",
    "            if result[i][j] == 1:\n",
    "                feature_in = pd.Series(result[i])\n",
    "\n",
    "                feature_out = result[i].copy()\n",
    "                feature_out[j] = 0\n",
    "                feature_out = pd.Series(feature_out)\n",
    "\n",
    "            else:\n",
    "                feature_out = pd.Series(result[i])\n",
    "\n",
    "                feature_in = result[i].copy()\n",
    "                feature_in[j] = 1\n",
    "                feature_in = pd.Series(feature_in)\n",
    "            \n",
    "            # -----------------------------------------------------------------\n",
    "            feature_in = feature_in[feature_in == 1].index.tolist()\n",
    "            feature_out = feature_out[feature_out == 1].index.tolist()\n",
    "            \n",
    "            if len(feature_out) > 0:\n",
    "                try:\n",
    "                # 实际使用中可能会因为进入模型的变量过多产生奇异矩阵，因此使用try来过滤这种特殊情况\n",
    "                    bic_out = temp * ic(data.loc[:,feature_out],criterion)\n",
    "                    try:\n",
    "                        bic_in = temp * ic(data.loc[:,feature_in],criterion)\n",
    "                        \n",
    "                        p = 1 / (1 + np.exp(bic_in - bic_out))\n",
    "                        \n",
    "                        rand = random.random()\n",
    "                        if rand < p:\n",
    "                            result[i][j] = 1\n",
    "                            criterion_value[i][j] = bic_in / temp\n",
    "                            \n",
    "                        else:\n",
    "                            result[i][j] = 0\n",
    "                            criterion_value[i][j] = bic_out / temp\n",
    "                            \n",
    "                            record[i][j] = 1\n",
    "                        print(str(100 * i / turns) + '%' + '_succeed')\n",
    "                    except:\n",
    "                        criterion_value[i][j] = bic_out / temp\n",
    "                        result[i][j] = 0\n",
    "                        record[i][j] = 0\n",
    "                        print(str(100 * i / turns) + '%' + '_failed')\n",
    "\n",
    "                except:\n",
    "                    try:\n",
    "                        bic_in = temp * ic(data.loc[:,feature_in],criterion)\n",
    "                        criterion_value[i][j] = bic_in / temp\n",
    "                        result[i][j] = 1\n",
    "                        record[i][j] = 1\n",
    "                    except:\n",
    "                        result[i][j] = 0\n",
    "                        record[i][j] = 0\n",
    "                        continue\n",
    "                    print(str(100 * i / turns) + '%' + '_failed')\n",
    "                    \n",
    "            else:\n",
    "                try:\n",
    "                    bic_in = temp * ic(data.loc[:,feature_in],criterion)\n",
    "                    criterion_value[i][j] = bic_in / temp\n",
    "                    result[i][j] = 1\n",
    "                    record[i][j] = 1\n",
    "                    print(str(100 * i / turns) + '%' + '_succeed')\n",
    "                except:\n",
    "                    result[i][j] = 0\n",
    "                    record[i][j] = 0\n",
    "                    print(str(100 * i / turns) + '%' + '_failed')\n",
    "                \n",
    "            k += 1\n",
    "        print(str(100 * i / turns) + '%')\n",
    "        \n",
    "    result = pd.DataFrame(result)\n",
    "    criterion_value = pd.DataFrame(criterion_value)\n",
    "    record = pd.DataFrame(record)\n",
    "    confidence = result.iloc[:,1:].sum(axis = 1) / turns\n",
    "    return [confidence, result, criterion_value, record]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3eec9cac",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = pd.read_csv(\"data.csv\")\n",
    "a.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "464a37ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = a.sample(n=20000, frac=None, replace=False, weights=None, random_state=0, axis=0)\n",
    "data = df.iloc[:,2:]\n",
    "data.index = list(range(data.shape[0]))\n",
    "sample = df.iloc[:,1]\n",
    "sample.index = list(range(sample.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18a40c19",
   "metadata": {},
   "outputs": [],
   "source": [
    "turns = 400 # 模拟轮数\n",
    "temp = 0.5 # 温度系数\n",
    "criterion = 'bic' # 信息准则，可以是aic 或者 bic\n",
    "confidence, result, criterion_value, record = gibbs_sampling(data,sample,turns,temp,criterion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfa4ae1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "result.to_csv(\"result.csv\",index=True)\n",
    "confidence.to_csv(\"confidence.csv\",index=True)\n",
    "record.to_csv(\"record.csv\",index=True)\n",
    "criterion_value.to_csv(\"criterion_value.csv\",index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28dd66e5",
   "metadata": {},
   "source": [
    "# 特征筛选"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdab43c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "turns = 400\n",
    "confidence = result.iloc[:,1:].sum(axis = 1) / turns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7073b1ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = {}\n",
    "for i in np.linspace(0, 1):\n",
    "    feature = list(confidence[confidence >= i].index)\n",
    "    for j in range(len(feature)):\n",
    "        feature[j] = feature[j][:-1]\n",
    "    feature = list(set(feature))\n",
    "    X = df.loc[:, feature]\n",
    "    X_train, X_test, y_train, y_test = train_test_split(X,\n",
    "                                                        y,\n",
    "                                                        random_state=0,\n",
    "                                                        test_size=0.2)\n",
    "    model = LogisticRegression()\n",
    "    model.fit(X_train, y_train)\n",
    "    y_proba_dtc_logi = model.predict_proba(X_test)\n",
    "    test_preds = model.predict(X_test)\n",
    "    auc_dtc = roc_auc_score(y_test, y_proba_dtc_logi[:, 1])\n",
    "    a[i] = auc_dtc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7bb81f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "b = zip(a.values(),a.keys())\n",
    "s = sorted(b)\n",
    "s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10934c78",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_list = []\n",
    "index = []\n",
    "for i in np.linspace(0,1):\n",
    "    number = len(confidence[confidence>=i])\n",
    "    index.append(i)\n",
    "    num_list.append(number)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "340a08da",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b536f589",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = []\n",
    "y = []\n",
    "for i in range(len(s)):\n",
    "    y.append(s[i][0])\n",
    "    x.append(s[i][1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fac130ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(7,5))\n",
    "plt.plot(x,y,color='b')\n",
    "plt.yticks(fontproperties = 'Times New Roman', size = 20)\n",
    "plt.xticks(fontproperties = 'Times New Roman', size = 20)\n",
    "plt.ylabel('AUC', fontdict={'family' : 'Times New Roman', 'size'   : 15})\n",
    "plt.xlabel('confidence', fontdict={'family' : 'Times New Roman', 'size'   : 15})\n",
    "plt.savefig(\"AUC遍历.png\",dpi = 300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5450d066",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(5,5))\n",
    "plt.plot(index,num_list,label='特征数量',color='b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "759e4237",
   "metadata": {},
   "outputs": [],
   "source": [
    "feature = list(confidence[confidence>=0.7].index)\n",
    "for i in range(len(feature)):\n",
    "    feature[i] = feature[i][:-1]\n",
    "feature = list(set(feature))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ee0b5ee",
   "metadata": {},
   "source": [
    "# 评价指标对比"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3d60570",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_excel(\"18WOE.xlsx\",header=0)\n",
    "data.to_csv(\"18WOE.csv\",index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a58a3d2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_before = data.iloc[:,:-1]\n",
    "Y_before = data.iloc[:,-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb07ff20",
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import Counter\n",
    "print(Counter(Y_before))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "502052a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from imblearn.over_sampling import SMOTE\n",
    "smo = SMOTE(sampling_strategy={0: 261784,1:65446},random_state=0)\n",
    "X_smo, y_smo = smo.fit_resample(X_before, Y_before)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e76e6d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_smo.to_csv(\"logstic_smo_x.csv\",index=True)\n",
    "y_smo.to_csv(\"logstic_smo_y.csv\",index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58f577dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_before, X_test_before, y_train_before, y_test_before = train_test_split(X_smo, y_smo, random_state = 30, test_size = 0.3)\n",
    "\n",
    "unique_, count_ = np.unique(y_train_before, return_counts=True)\n",
    "for unique, count in zip(unique_, count_):\n",
    "    print(\"训练集中类别为'{0}'的样本有{1}个\".format(unique, count))\n",
    "\n",
    "print()\n",
    "unique_, count_ = np.unique(y_test_before, return_counts=True)\n",
    "for unique, count in zip(unique_, count_):\n",
    "    print(\"测试集中类别为'{0}'的样本有{1}个\".format(unique, count))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77995582",
   "metadata": {},
   "outputs": [],
   "source": [
    "model=LogisticRegression(C=2, solver='saga', class_weight={1: 1, 0: 3})\n",
    "model.fit(X_train_before, y_train_before)\n",
    "y_proba_before=model.predict_proba(X_test_before)\n",
    "test_preds = model.predict(X_test_before)\n",
    "acc_dtc = accuracy_score(y_test_before,test_preds)\n",
    "kappa_dtc = cohen_kappa_score(y_test_before,test_preds)\n",
    "auc_dtc = roc_auc_score(y_test_before, y_proba_before[:,1])\n",
    "fpr_dtc, tpr_dtc, thresholds_roc_dtc = roc_curve(y_test_before,y_proba_before[:,1])\n",
    "ks_dtc = max([trp_ - fpr_ for trp_, fpr_ in zip(tpr_dtc, fpr_dtc)])\n",
    "results = pd.DataFrame()\n",
    "results['准确率'] = [acc_dtc]\n",
    "results['AUC'] = [auc_dtc]\n",
    "results['K-S值'] = [ks_dtc]\n",
    "results['kappa'] = [kappa_dtc]\n",
    "results.index = ['dtc']\n",
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7898f9f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_excel(\"WOE30.xlsx\")\n",
    "data.to_csv(\"WOE30.csv\",index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "544df738",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = data.iloc[:,:-1]\n",
    "Y = data.iloc[:,-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fa378f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import Counter\n",
    "print(Counter(Y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2641d0ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "from imblearn.over_sampling import SMOTE\n",
    "smo = SMOTE(sampling_strategy={0: 261784,1:65446},random_state=0)\n",
    "X, y = smo.fit_resample(X, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "271f96f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "X.to_csv(\"logstic_smo_x30.csv\",index=True)\n",
    "y.to_csv(\"logstic_smo_y30.csv\",index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56afb3b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = X\n",
    "X = df.loc[:,feature]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e070a708",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X,y,random_state = 0,test_size = 0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02621117",
   "metadata": {},
   "outputs": [],
   "source": [
    "model=LogisticRegression(C=1, solver='saga', class_weight={1: 1, 0: 1})\n",
    "model.fit(X_train, y_train)\n",
    "y_proba_dtc_logi=model.predict_proba(X_test)\n",
    "test_preds = model.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1462806",
   "metadata": {},
   "outputs": [],
   "source": [
    "acc_dtc = accuracy_score(y_test,test_preds)\n",
    "kappa_dtc = cohen_kappa_score(y_test,test_preds)\n",
    "auc_dtc = roc_auc_score(y_test, y_proba_dtc_logi[:,1])\n",
    "fpr_dtc, tpr_dtc, thresholds_roc_dtc = roc_curve(y_test,y_proba_dtc_logi[:,1])\n",
    "ks_dtc = max([trp_ - fpr_ for trp_, fpr_ in zip(tpr_dtc, fpr_dtc)])\n",
    "results = pd.DataFrame()\n",
    "results['准确率'] = [acc_dtc]\n",
    "results['AUC'] = [auc_dtc]\n",
    "results['K-S值'] = [ks_dtc]\n",
    "results['kappa'] = [kappa_dtc]\n",
    "results.index = ['dtc']\n",
    "results"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": false,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
